Signatures of unconventional pairing in near- vortex electronic structure of LiFeAs 
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A major question in Fe-based superconductors remains the structure of the pairing, in particular 
whether it is of unconventional nature. The electronic structure near vortices can serve as a plat- 
form for phase-sensitive measurements to answer this question. By solving Bogoliubov-de Gennes 
equations for LiFeAs, we calculate the energy-dependent local electronic structure near a vortex for 
different nodeless gap-structure possibilities. At low energies, the local density of states (LDOS) 
around a vortex is determined by the normal-state electronic structure. However, at energies closer 
to the gap value, the LDOS can distinguish an anisotropic from a conventional isotropic s-wave gap. 
We show within our self-consistent calculation that in addition, the local gap profile differs between a 
conventional and an unconventional pairing. We explain this through admixing of a secondary order 
parameter within Ginzburg-Landau theory. In-held scanning tunneling spectroscopy near vortices 
can therefore be used as a real-space probe of the gap structure. 

PACS numbers: . . . 



I. INTRODUCTION 

The gap structure in the iron-based superconductors 
and its possible unconventional nature is still a key issue 
in the field four years after their discovery. In most com- 
pounds of the family, the pairing is believed to be of the 
so-called 8 ± type, for which the order parameter changes 
sign between the electron-like and the hole-like Fermi 
surfaces. 1 However, a major difficulty in distinguishing 
such an unconventional pairing state from a trivial s- 
wave gap is that both states are nodeless and transform 
trivially under all the symmetry operations of the mate- 
rial's point group. As the experimental probes that are 
usually used to distinguish various gap structures, such 
as phase- sensitive probes, are not Fermi pocket specific, 
an unambiguous evidence of the unconventional s ± pair- 
ing remains evasive. 

One route to accessing phase information using a phase 
insensitive probe would be through vortex bound states, 
as a vortex introduces a spatial texture to the super- 
conducting order parameter. Advancements in in- field 
scanning tunneling spectroscopy (STS) enables study of 
vortex bound states. Indeed, a recent scanning tunnel- 
ing spectroscopy experiment on LiFeAs under magnetic 
field has shown intriguing energy dependence in the spa- 
tial distribution of local density of states (LDOS) near a 
vortex. 2 At zero bias, the LDOS shows a four-fold star 
shape with high-intensity 'rays' along the Fe-As direc- 
tion, and with increasing energy it exhibits 'splitting of 
the rays'. The most striking features emerge at energies 
E w A/2 as hot spots, at the intersections of two rays. 

Motivated by these observations, we present a study of 
the electronic structure and signature of unconventional 
pairing near a vortex, within the Bogoliubov-de Gennes 
(BdG) framework. By imposing a gap structure and solv- 
ing the BdG Hamiltonian, we show that the isotropic s- 
wave and s ± -wave pairing result in different spatial dis- 
tributions of the LDOS at energies approaching the gap 
value. In particular, we find s ± -wave pairing to yield 



observed hot spots. 

For a more direct evidence of unconventional s^-wave 
pairing, we propose detecting spatial distribution of the 
gap around a vortex, based on our self-consistent solu- 
tions to the BdG equations. In a conventional isotropic 
s- wave superconductor, a vortex introduces a singular 
point in space around which the phase of the order 
parameter winds. In addition, the gap amplitude is 
completely suppressed at the core of the vortex, but 
smoothly and isotropically recovers away from the core. 
However, in an unconventional superconductor, a vor- 
tex can also induce a secondary order parameter in its 
vicinity. 3 9 . This is possible since a symmetry distinc- 
tion between, for instance, a gap with s ± ~ cos k x cos k y 
and d xy ~ sin k x sin k y structure is only meaningful in 
a translationally invariant system. Due to the induced 
secondary order parameter (which diminishes far from 
the vortex), the gap recovery should show a strong angu- 
lar dependence. Detection of such anisotropy will be an 
unambiguous evidence of unconventional pairing. 

The rest of the paper consists of the following: In sec- 
tion II, we present our microscopic model. In section III, 
we describe our methods. In section IV, we present the 
results of BdG calculations and insight from correspond- 
ing Ginzburg-Landau theory. In section V, we summarize 
the our findings and remark on future directions. 

II. MODEL 

We describe LiFeAs in the superconducting state with 
the (mean-field) BdG Hamiltonian 

^ BdG = E*!( A |' 

Here, = (c^,cj^) T is a Nambu spinor, and q s (cj s ) 
annihilates (creates) an electron at lattice site i with spin 
s. The kinetic part Uj denotes hopping between sites i 
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FIG. 1. (a) Fermi surface of LiFeAs with three hole- like 
pockets around the Y point and the two electron-like pockets 
around the M point. In this work, we focus on the 7 band 
whose Fermi surface is shown as a bold line, (b) Sketch of the 
three gap functions with s-, s ± -, and d^-wave momentum 
structure around the 7-band Fermi surface. 



and j, and A^- are the (bond) gap functions. To solve 
H BdG self-consistently, the gap functions have to satisfy 
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(2) 



where Vij < is the attractive interaction between sites 
i and j in the singlet channel and (•) denotes the thermal 
expectation value. Restricting the interaction Vij to a 
specific form constrains the momentum structure of the 
gap function, since A^- ^ only if Vij ^ 0. In the 
uniform case, on-site attraction Vij = USij leads to spin- 
singlet s-wave gap 



A.(k) =A°, (3a) 
while next-nearest-neighbor (NNN) attraction Vij = 



V'S^j^ allows for the singlet gap functions of s 
A s ± (k) =4A^± cos k x cos k y 



± 



and d xy form 



A d xy (k) =4A° dxy sin k x sin k y 



(3b) 



(3c) 



in the reciprocal space. 

In this paper, we focus on a single band: the so-called 
7 band. The Fermi surface of the 7 band is shown in 
Fig. 1(a) (solid line) along with four other bands (dashed 
lines) crossing the Fermi level in LiFeAs. 10 13 This choice 
is motivated by the 7 band's smaller gap value. 11 ' 14 More- 
over, this band mainly stems from the (in-plane) d xy 
orbitals, and hence shows little k z dependence. 12 Our 
choice of tij, guided by the experimental observations on 
the 7 pocket, 11,14 ' 15 is t = — 0.25eV for nearest-neighbor 
hopping, t' = 0.082eV for next-nearest-neighbor hopping, 
and ta = n = 0.57eV for the chemical potential. No- 
tice the highly anisotropic (square-like) Fermi surface of 



this band shown in Fig. 1(a). At low energies, these flat 
(quasi-one-dimensional) parts of the electronic structure 
dominate over the small anisotropy of the «s ± gap [see 
Fig. 1(b)] in determining the near- vortex LDOS. 

We introduce vortices into the system through two ap- 
proaches. First, we fix the non-zero gap functions for 
both on-site 8-wave pairing and NNN s^-wave pairing 
by hand to 



=A°tanh(|r^|/0e il 



(4) 



where the vector r^- points to the midpoint of sites i 
and j, and Oij is its azimuthal angle measured from the 
Fe-Fe direction. This therefore corresponds to a single 
vortex suppressing the gap at the core located at the 
origin and with phase winding around it. Second, we also 
perform self-consistent calculations, where we induce the 
formation of vortices by applying a magnetic field Hz 
on the system. Assuming minimal coupling between an 
electron and the field, the hopping between sites i and j 
acquires a Peierls phase 
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A(r) • dr, 



(5) 



(6) 



where §0 = h/2eis the magnetic fluxoid and is the vec- 
tor pointing to site i. We assume uniform magnetic field 
H and write the vector potential in the Landau gauge 
A(r) = -Hy±. 

From the self-consistent solution A^-, we can define 
local gap order parameters of different symmetries. For 
an on-site interaction, the local s-wave order parameter 
is defined as 



A s (r) =A r 



(7a) 



Note that from here on, we use r without any site index to 
denote both a lattice site and the corresponding vector in 
units of the lattice constant a^. With NNN interaction, 
we define local order parameters of s ± 

A s ±(r) = [A r+ (i ?1 ) 5r + A r+ ( l5 _i) 5r 

+ A r+( _ l5 _ 1)5r + A r+( _ M) , r ]/4 (7b) 

and d xy form 

^d xy (r) = [A r+ (i 5 i) ?r - A r+ (i _i), r 

+ A r+( _ l5 _ 1))r - A r+( _ M) , r ]/4, (7c) 

where 

A rr , =A rr ,e-^ r ' r '\ (8) 

ensures that order parameters of different symmetries do 
not mix under magnetic translations. Note that for the 
uniform case, A s (r) = A°, A s ±(r) = A° ± and (r) = 
A° ofEqs. (3a)-(3c). 
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III. METHOD 

In this section, we elaborate on our two approaches to 
solve the BdG equations and obtain the local density of 
states near a vortex. For both, diagonalizing the Hamil- 
tonian H BdG in Eq. (1) for a system of size (N x ,N y ) is 
computationally the most expensive part. 

A. Non-Self-Consistent Approach 

For the non-self-consistent calculation, we impose a 
gap function in the form given by Eq. (4) and find the 
low lying eigenvalues and eigenstates of l~L BdG using the 
Lanczos algorithm. We suppress the low energy states 
from forming at the boundary by applying a large poten- 
tial to the boundary sites. The LDOS can be calculated 
from the eigenvalues E n and eigenstates [u n (r) , v n (r)] as 

N(r, E)ocY^ \u n (v)\ 2 8(E - E n ) + \v n (r)\ 2 5(E + E n ). 

n 

(9) 

Since we are not interested in the absolute value of the 
LDOS but rather in the spatial profile at a given energy, 
we normalize the LDOS such that for a given energy E, 
the maximum value of iV(r, E) is unity. 

B. Self-Consistent Approach 

For the self-consistent calculation, we assume some ini- 
tial gap functions and use the eigenvalues and eigenvec- 
tors of Eq. (1) to calculate the gap functions given by 
Eq. (2). We proceed iteratively until self-consistency is 
achieved. In diagonalizing H BdG , we can no longer make 
use of the crystal momentum basis to simplify the prob- 
lem since the Peierls phase factor prevents the kinetic 
part of the Hamiltonian from commuting with the ordi- 
nary lattice translation operator Tr. However, the ki- 
netic part commutes with the magnetic translation oper- 
ator 

f R = e-^ A(R)r T R . (10) 

for a magnetic lattice vector R whose unit cell contains 
two magnetic fluxoids. 

The pairing term in general does not commute with 
Tr. Nevertheless, when vortices form a lattice, Tr com- 
mutes with the pairing term if R is a vector of a vor- 
tex sublattice containing every other vortex. Since we 
focus on the electronic structure near a single vortex, 
we expect the shape of the vortex lattice to have little 
effect on our results. Therefore, we make an arbitrary 
choice for its primitive vectors to be T x x and L y y, such 
that R forms a rectangular lattice R = (m x L x ,m y L y ), 
where m a = • • • M a - 1 and M a N a /L a . 16 Note 
that periodic boundary conditions in the Landau gauge 
A(r) = —Hy-k requires the total magnetic flux through 



the system to be an integer multiple of 2<fr N x . In addi- 
tion, one magnetic unit cell contains a magnetic flux of 
2<3>ch i- e - H = 2<I>o/ L x L y . We satisfy these two require- 
ments by choosing M x = L y , M y = L x . 
Working with the magnetic Bloch states 

* k (r) = ^e- ikR f R *(r)f^ 1 (11) 

R 

allows us to block diagonalize the Hamiltonian 

nBdG = ww £ £ r '^^- ( 12 ) 

x y k r,r' 

The indices k and r are from here on defined in the mag- 
netic Brillouin zone and magnetic unit cell, respectively, 
i.e., 

r=(£ x ,£ y ) £ a = 0---L a -l. (13b) 

By diagonalizing the matrices of dimension 2L x L y x 
2L x L y , we can compute the eigenstates and eigenenergies 
of T~L BdG . These are then used to calculate A^- with 
Eq. (2) closing the self-consistency cycle. 

We use the self-consistent solution to calculate the 
local order parameters of specific symmetries as defined 
in Eqs. (7a)- (7c) and also the LDOS of the electronic 
degrees of freedom, as defined in Eq. (9). 

IV. RESULTS 
A. Non-Self-Consistent Approach 

Figure 2 shows the near-vortex LDOS calculated by 
diagonalizing H BdG in Eq. (1) with fixed gap functions 
as given by Eq. (4) on a system of dimension (N X: N y ) = 
(301,301). We choose realistic values of the parameters 
for the coherence length £ = 16. 4ao 17,18 and gap values 
A® = 3meV for on-site pairing and A° ± = 1.5meV for 
NNN pairing. 11 ' 14 ' 15 

We can interpret the vortex bound states in this non- 
self-consistent BdG calculation as bound states in a trap 
potential given by Eq. (4), where only states around the 
normal state Fermi surface constitute the bound states. 
There are then two sources of anisotropy: anisotropic, 
quasi-one-dimensional low-energy properties of the nor- 
mal state, and an anisotropic gap, both defined in the 
momentum space in the translationally-invariant system. 

We first focus on the low energy LDOS [Figs. 2(a) and 
(b)], where the normal-state properties dominate. For 
both types of pairing, the bound state is located at the 
center of the trap. Since the Bloch states making up this 
bound state have two main velocities due to the quasi- 
one-dimensional parts of the Fermi surface, the bound 
state only extends in these two directions out of the trap, 



4 





FIG. 2. Local density of states near a vortex for non-self-consistent calculation with the gap function given by Eq. (4). 
The value iV(r, E) has been normalized such that the maximum value in each map is unity, (a) is the LDOS at the lowest 
bound state energy with on-site pairing with A = 3meV, and (d) is at higher energy, (b) and (e) are with NNN pairing with 
A = 1.5meV. (c) and (f) are the near-vortex LDOS maps observed by Hanaguri et al 2 . 



(b) 



resulting in the rays in Fig. 2(a) and (b). The gap is 
suppressed near the vortex center, and its anisotropy is 
of little importance. 

At higher energies, on the other hand, the bound state 
is located away from the vortex core. The quasi-one- 
dimensionality of the Fermi surface allows for localiza- 
tion in one direction and extended parts in the other 
direction. This leads to a square-like inner ring in the 
LDOS for both cases [Figs. 2(c) and (d)]. The difference, 
however, results from the anisotropy of the gap function. 
While the isotropic s-wave gap is analogous to a potential 
that is independent of the momentum, the anisotropic 
gap is a trap for which different states around the Fermi 
surface experience different potentials depending on their 
momenta. With the gap function of s ± form, the quasi- 
one-dimensional portion of the Fermi surface experiences 
a stronger trap potential, leading to a suppression of its 
contribution to the bound-state wave function. As a re- 
sult, the bound state exhibits pronounced isolated seg- 
ments, 'hot spots', within the inner ring that point in 
the Fe-Fe direction, as shown in Fig. 2(d). 



B. Self-Consistent Approach 
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FIG. 3. Local density of states near a vortex for self- consistent 
calculation. Here again, N(r, E) is normalized, (a) is the 
LDOS at the lowest bound state energy with on-site attraction 
U — — eV, and (c) is at higher energy, (b) and (d) are with 
NNN attraction V' = -0.3eV. 



Figure 3 shows the results from the self-consistent cal- 
culation. We compare two pairing interactions - on-site 
attraction U = -0.35eV, and NNN attraction V = 
— 0.3eV - for a system with magnetic unit cell of dimen- 
sions (L x ,L y ) = (19,38). This corresponds to a full lat- 
tice size of (N x ,N y ) = (38 x 19,19 x 38). In zero field, 



the two cases lead to a uniform superconducting gap of 
A® = 27meV and A^ ± = lOmeV, respectively. We have 
chosen U and V such that the coherence length £ oc A -1 
is small compared to the inter-vortex spacing. This al- 
lows us to focus on a nearly isolated vortex within the 
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FIG. 4. Spatial distribution of different symmetry compo- 
nents of order parameters, (a) A s (r) for on-site attraction 
U = -0.35eV. (b) A s ±(r) and (c) A dxy (r) for NNN pairing 
V' = — 0.3eV. The values have been normalized by the value 
of the dominant order parameter in the absence of magnetic 
field for each case: A° for (a), and A° s ± for (b), (c). The 
contours in red indicate points at which the amplitude of the 
order parameter is 0.825 for the innermost, and 0.925 for the 
outermost contours, with equal intervals between the con- 
tours in between. Note that the color-scale for Ad (r) much 
smaller than for A s (r). 



computationally feasible size of the magnetic unit cell. 
Although the resulting gap values are an order of mag- 
nitude larger than what is known experimentally, this 
should not affect the validity of the results in a qualita- 
tive manner. Both at low energy and at higher energy 
close to the gap value, we observe features that quali- 
tatively agree with the results obtained in the previous 
section. 

The self-consistent calculation also allows us to study 
the local order parameters of a given structure near a 
vortex. Unlike for the on-site attraction, where A s (r) is 
the only allowed gap function, order parameters of dif- 
ferent symmetries can mix near a vortex for NNN at- 
traction. A near- vortex map of A s (r) for on-site pairing 
shown in Fig. 4(a) indeed shows isotropic healing of the 
order parameter away from the vortex core. On the other 
hand, for the NNN attraction which leads to uniform s- 
wave pairing in zero-field, the secondary order parameter 
Ad xy (r) is induced near the vortex. Coupling between 
this secondary order parameter and the primary A s ± (r) 
leads to strong angular variation of both components as 
can be seen in Figs. 4(b) and (c). 

To gain further insight into the admixing of a sec- 
ondary order parameter near a vortex for the anisotropic 
pairing, we analyze the Ginzburg-Landau free-energy 
density. While an on-site attraction only allows for an 



isotropic s-wave component, the NNN pairing allows for 
both s- and d-wave components. Hence, the symmetry- 
based free-energy density corresponding to the NNN- 
anisotropic-pairing case is 

/ =a s \s\ 2 + a d \d\ 2 + AM 4 + W + fo\s\ 2 \d\ 2 
+ f3 4 (s* 2 d 2 + c.c.) + 7s \Ds\ 2 + ld \Dd\ 2 
+ j u (D x sD y d* + D y sD x d* + c.c), (14) 

where s and d are shorthands for s(r) and d(r), the or- 
der parameter fields for the s ± and d xy gaps, respec- 
tively, and Di = di — ieAi is the covariant derivative. 
The fields s(r) and d(r) can be thought of as A s / S ±(r) 
and Ad xy (r) after coarse graining. This type of admix- 
ing near a vortex has previously been studied in the con- 
text of cuprates. 3 919 As the large halo around vortices 
in cuprates 20 hindered the observation of this admixing, 
LiFeAs presents an opportunity for this observation. 

The spatial variation of the secondary component d xy 
in Fig. 4(c) is largely due to the derivative coupling, the 
term proportional to 7^ in Eq. (14). For \s\ ^> \d\ and 
I Ds I ^> \Dd\ the spatial structure of the d xy component is 
determined largely by the structure of the s-wave com- 
ponent. Minimizing Eq. (14) with respect to d(r) and 
keeping only terms up to linear order in d(r), we find 

-j d D 2 d + a d d + f3 3 \s\ 2 d + pAS 2 d* =j u (D x D y + D y D x )s. 

(15) 

Hence, the curvature in the leading s- wave component 
will induce the secondary (d xy ) component. Now con- 
sider a single isolated vortex. As s(r) is recovered at 
the length scale of the coherence length £ away from the 
core of the vortex, we expect a large d(r) due to coupling 
to the large curvature of s(r) at this distance. Since 
£ = Hvf/ttA ~ 3.0ao for the uniform gap value with 
V = — 0.3eV, this is in agreement with the positions of 
the maxima of d(r) in Fig. 4(c) as a function of |r| set- 
ting the vortex core at the origin. We can also explain the 
angular variation and the form of the anisotropy of d(r) 
in this framework. If we assume s(r) = f(r)e ze with a 
slowly changing f(r) and the azimuthal angle measured 
from the Fe-Fe direction, we find from Eq. (15) 

d(r) - d x d y s(r) - e~ ie {l + 3e 4 ^), (16) 

ignoring the phase due to the magnetic field. The struc- 
ture of the derivative hence gives rise to a four-fold 
anisotropy, which explains the fact that \d(r)\ is maxi- 
mum in the Fe-Fe direction, while it is suppressed along 
the 45° direction. Coupling to d(r) gives then in turn 
cause for the four-fold anisotropy in s(r). 

V. CONCLUSION 

We have contrasted effects of anisotropic s^-wave 
(NNN pairing) and isotropic s-wave (on-site pairing) on 
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the near- vortex local density of states in LiFeAs, by 
solving Bogoliubov-de Gennes equations both non-self- 
consistently and self-consistently. We find qualitative 
changes in the geometric distribution of the density of 
states as a function of energy. At low energies, the 
anisotropy of the vortex bound state, and hence the 
LDOS, are determined by the normal-state low energy 
electronic structure, independently of the gap structure. 
The different pairing structures, however, lead to qual- 
itatively different LDOS distribution at higher energies: 
While the isotropic s-wave shows a square-like feature of 
roughly equal intensity, four 'hot spots' develop in the 
case of an (anisotropic) ,s ± -wave gap. Indeed, our re- 
sults for the latter case qualitatively agree with recent 
experiments by Hanaguri et al. 2 . 

From a self-consistent treatment we further find a dif- 
ference in the recovery of the order parameter away from 
the vortex core: a pronounced angular dependence of 
the s^-wave gap compared to isotropic behavior for the 
s- wave gap. Employing a Ginzburg-Landau analysis, we 
explain this difference through admixing of a secondary 
order parameter supported by the NNN interaction. De- 
tection of the anisotropy or even the secondary order pa- 
rameter would be a strong proof of the unconventional 
nature of the pairing. 

Here we focused on the 7 band with interest in low 
energy properties, as this is the band with the smallest 
gap. 11,14 Hence for features at energies much less than the 
gap energy at the electron pockets, we expect our calcu- 
lation to capture salient features of in-field STS exper- 
iments. The comparison between the calculated LDOS 



and the results in Ref. 2 confirms this conjecture. 

In closing we note that our calculation captures 
Friedel-like oscillations, frequently referred to as quasi- 
particle interference (QPI), due to vortices. QPI in the 
presence of vortices was successfully used to access phase 
information with STS in cuprates. 21 Recent in-field QPI 
experiments on FeSe have been interpreted to be con- 
sistent with an s ± scenario when a vortex is treated as 
a magnetic scatterer for BdG quasiparticles. 22 However, 
a vortex is at once a point of gap suppression, a point 
with magnetic flux and a point around which the order- 
parameter phase winds. While we treated vortices faith- 
fully in this work, we could not investigate effects of inter- 
pocket sign change as we only considered one pocket. An 
extension of the present work with the full band struc- 
ture would be necessary to work out what to expect for 
different order-parameter possibilities, especially how the 
phase difference between different pockets affects in-field 
QPI. 
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